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Abstract 

The Boeder differential equation is solved in this work over a wide range of a, yielding the prob- 
ability density functions (PDF), that describe the average orientations of rod-like macromolecules 
in a flowing liquid. The quantity a is the ratio of the hydrodynamic shear rate to the rotational 
diffusion coefficient. It characterises the coupling of the motion of the macromolecules in the 
hydrodynamic flow to their thermal diffusion. Previous analytical work is limited to approximate 
solutions for small values of a. Special analytical as well as numerical methods are developed in 
the present work in order to calculate accurately the PDF for a range of a covering several orders 
of magnitude, 10~^ < a < 10^. The mathematical nature of the differential equation is revealed 
as a singular perturbation problem when a becomes large. Scaling results are obtained over the 
differential equation for a > 10^. Monte Carlo Brownian simulations are also constructed and 
shown to agree with the numerical solutions of the differential equation in the bulk of the flowing 
liquid, for an extensive range of a. This confirms the robustness of the developed analytical and 
numerical methods. 
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I. INTRODUCTION 



The most widely used experimental technique for the observation of rotational diffusion 
of macromolecules in hydrodynamic flow is that corresponding to birefringence [1,2]. The 
authors in [2] , for example, have observed the existence of a depletion layer in dilute xanthan 
solutions subjected to simple shear flow, using the technique of evanescent induced fluores- 
cence. The motion and average orientations of macromolecules at dilute concentrations in 
a flowing liquid are determined by two conflicting forces, the first is a hydrodynamic force 
stemming from shear flow, and the other thermal, originating from Brownian rotational dif- 
fusion. 

Boeder [3] is the first to have studied this problem in the bulk of a flowing liquid from the 
theoretical point of view. He suggested, without the use of the Langevin formalism or the 
Fokker-Planck equation [4,5], an ordinary differential equation, which governs a probability 
distribution function P{9) describing the average orientations of the macromolecules: 

P" + {asin'^{e)Py = (1) 

This differential equation is derived for the motion of macromolecular rod-like particles, 
of negligible cross-sectional area, in the plane of the flow, without any boundary conditions. 
The angle h describes the orientations of the macromolecules with respect to a reference 
direction. The quantity a is the ratio 



where 7 is the shear rate in the hydrodynamic flow, and Dbt is a diffusion coefficient 
governing the Brownian rotational motion of the rod. The differing nature of the macro- 
molecules and macromolecular polymers in their solutions leads to differing values of Dsr, 
whereas the hydrodynamic flow of the solution is characterised by given values of 7. As has 
been pointed out [3], a solution in closed form for this ordinary differential equation cannot 
be obtained. Boeder gave a series solution valid only for small values of a. Although some 
improvements have been made, on one hand to remove restrictions on the cross-sectional 
areas of the macromolecules, and on the other to be able to consider rotational diffusion in 
three dimensions [6,7], the available solutions are still limited to small values of a. To the 
best of the authors' knowledge there does not exist in the literature a general approach to 
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solve the Boeder differential equation for arbitrarily large values of a, or for a wide range 
of its values. It is the purpose of this work to provide such an approach, using different 
analytical and numerical methods. Another interest is to compare the analytical approach 
with Brownian simulations of the PDF, under the same conditions in the bulk of a flowing 
liquid. This comparison is useful because it can confirm the robustness of the analytical 
methods. It also provides a necessary limiting bulk condition for similar simulations near 
solid surfaces where one is often led to making reasonable but unverifiable assumptions con- 
cerning the touching collision of the macromolecule with the sohd surface, [4,5,8]. Although 
simulations have been developed in the past, as for polymers in plane bulk Poiseuille flow [8], 
and for the dispersion of rod-like particles in a bulk shear flow [9], there is an absence of a 
detailed comparison between the simulations and the analytical results in the bulk, precisely 
because exact solutions of this differential equation have not been available. In Section 2, 
we present an accurate analysis of the Boeder differential equation to obtain the probability 
distribution function (PDF), for a wide range of a. In Section 3, four numerical methods are 
proposed for the general solution of the differential equation in different intervals of a, and 
numerical results are presented to illustrate this. Scaling results are also presented in this 
section for large values of a. In Section 4, Monte Carlo Brownian simulations are presented 
for the PDF, and compared with the results of Section 3. The conclusions are given in 
Section 5. 2. 

II. ACCURATE ANALYSIS OF THE BOEDER EQUATION 

In this section, a procedure for the accurate numerical analysis of the Boeder ordinary 
differential equation and its associated probability distribution function, P{0), is given for a 
wide range of a. Turbulence effects are known to take place for values of a > 10^ . Although 
the Boeder differential equation ceases to apply in a strict physical sense beyond this limit, 
the present mathematical analysis is not limited by this, and numerical solutions may be 
calculated in our approach for values of a > 10^. P{9) is the solution of a second-order 
ordinary differential equation, Eq. (1). Since the Boeder PDF is periodic with a period 
equal to tt, as the macromolccular rods are indistinguishable when oriented at ^ or ^ -(- tt, 
the periodic boundary conditions 
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P(0) = P(7r) and P'(0) = P'{n) 
apply. In addition, the PDF has to be normahsed over the intervah 



(3) 



rp(e)de^i (4) 

Jo 

The determination of the Boeder PDF is consequently a constrained boundary value 
problem. Eq. (1), however, is not the only possible form for the Boeder equation, actually 
two other versions of the equation exist. The first takes the form of: 

P' + asin^(e)P (5) 

provided the initial value P{0) and the constant C are known. The second may be 
presented as: 



du/dO = V, (6) 
dv/dO = -asin{9)[sin{9)v + 2cos{9)u] (7) 

where u{e) = P{e),v{e) = P\9). 

Several technical procedures may be used to determine P{9) , depending on the different 
possible formulations of the problem, as an initial value problem or as a boundary value 
problem. In all cases two quantities, the constant C in Eq. (5), equal to P'(0) or v{Qi) , and 
the initial value P(0) or m(0), have to be evaluated from a for any of its given values. We 
have developed two methods, described below, to determine C and P(0). Firstly, a direct 
method based on the solution of Eq. (5), and secondly a minimisation method based on a 
multidimensional secant or Broyden method. The direct method to evaluate C and P(0) is 
as follows. The formal solution of Eq. (5) may be given generally as 

re 

P{9) = Cexp[a[sin{29)/2 - 9]/2] / exp[-a{sin{2x)/2 - x)/2]dx (8) 

J — CO 

The lower limit — oo is surprising since the problem is defined over the angular interval 
[0, tt] It is probably true to say, however, that the assumption to limit the space of stochastic 
events to this interval has retarded the analytical approach until now. We can analytically 
show that the lower limit — oo is the only possibility compatible with the boundary conditions 
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given by Eq. (3). Physically this implies that the stochastic space is cumulative over an 
infinite number of events and hence angles. Performing a change of variables, the solution 
may then be written as 



Pid) = (2C/a) expla sin(29)/A\ / exp(-x) expl a siniAxIa - 29)/A\dx (9) 

Jo 

The form in Eq. (9) is more stable numerically than the previous one of Eq. (8), since 
the numerically troublesome exp{6) and exp{—9) terms are avoided. Nevertheless when a 

increases the exp[a sin{9)/A] term will cause problems despite the bounded values of the 
sine function, forcing us to turn to other methods based on extrapolation techniques. The 
constant C is next determined from the normalisation condition of the Boeder PDF, as in 
Eq. (4), yielding 



C = a/ 



t/2 

whereas P{0) is given by 



/7r/2 roo 
d9 / dx exp{—x) exp[{a/2)sin{2x/a)cos{29 — 2x/a)] 
-7r/2 Jo 



(10) 



roo 

P(0)^(2C/a) exp(-x)exp[asin(4:x/a)/4\dx (11) 
Jo 

The PDF depends on a which we want to vary over several orders of magnitude. The 
difficulty in solving the differential equation stems from the fact that its nature may be 
modified when a increases, turning the problem into a singular perturbation one in a~^. 



III. NUMERICAL SOLUTIONS AND SCALING OF THE BOEDER EQUATION 

In order to cope with the wide range over which a may vary, we classify the various 
methods for solving the problem according to the value of a. Essentially, there are at least 
four methods to evaluate the PDF. 



Method a. For a in the range : [10 ^, 10^] . Direct calculation from the analytic solution 
of the first-order differential equation. 

Method b. For a in the range : [10^'^, 10^]. If C and P{0) are known, a simple ID 
Runge-Kutta method is used to solve the first-order differential equation. 
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Method c. For a. in the range : [10~^, 10^] . Find C and P(0) by a minimisation 
method based on a multidimensional secant or Broyden method, and do a 2D Runge-Kutta 

method to solve the system of two first-order differential equations, Eqs. (6) and (7). The 
condition min|P(0) — -P(7r)| comes from periodicity, and that of min| Jq P{6) — 1| is from 
the normalisation of the PDF. 

Method d. For For a in the range : [10^, 10^] and above. Calculate C and P(0) by an 
extrapolation method and solve Eq. (5), or Eqs. (6) and (7), with singular perturbation 
integration methods. To illustrate the different numerical methods used in our approach to 
solve the Boeder differential equation, and to calculate the PDF, some numerical results are 
presented, for a wide range of values of a. In Figs. 1-4, the continuous curves depict these 
PDF functions, for relatively high values of a —100, 1000, 2000, and 5000, respectively. The 
PDF results are normalised with respect to unity. The PDF for other values of a have been 
numerically calculated in a wide range, though they are not given here. The scaling results 
are presented next for large values of a. Eq. (10) can be bounded by replacing the cosine 
term in the exponential integrand by 1, and reducing the double integration appearing in 
the C denominator into a simpler one, namely: 

roo 

Dc — I expi—x)exp\:^ial2)sini2xla)\dx (12) 
Jo 

The sine term in Eq. (12) may be expanded for large values of a to yield Dc ~ ct^^^ 
for the upper bound case (-1- sign case), and Dc ~ of* for the lower bound case ( - sign 
case). Since C — a/2Dc ; we obtain the limits on the behaviour of C, namely C ~ a^/^ for 
the upper bound, case (-1- sign case), and C ~ for the lower bound case ( - sign case). 
Consider that C ~ a^/^ is an intermediate case in the space of stochastic events, when a is 
large. It is then possible to show that: 

P(0) ~ a^l^ 

P{Oma.) - (13) 

The proof is as follows. Expanding the argument of the exponential integrand in Eq. 
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(11), yields 



POD 

P(0) = (2C/a) / exp{-8x^/3a'^)dx (14) 
Jo 

The integral in Eq. (14) behaves as ~ a^^^, hence the result expressed in first part of Eq. 
(13). In order to prove the result of second part of Eq. (13), consider Eq. (5) at ^ = Omax, 
which may be written in this case as 

a sin^{emax)P{Omax) = C (15) 

For large values of a, the abscissa O^ax of the maximum of the Boeder PDF is relatively 
small, and Eq. (15) may be approximated to yield 

« OlaxPiOmax) = C (16) 

Since 9max is small, we can expand P{9max) around ^ = 0; such that 

P{Omax) = P{0) + 9l,,P" {9max)/2 + ... (17) 

P{9max) — -P(O) is hence accurate to second-order in 9max- This may be seen graphically, 
for example, in Figs. 2 and 3 for values of a— 1000 and 2000, and also for higher values. The 
result expressed in last part of Eq. (13) follows accordingly. Another direct consequence 
follows from Eq. (17) since it is now possible to write that 

9maxP{0) ~ 9maxP{9max) ~ COnstaut (18) 

which has a simple geometric interpretation. From the normalisation of the Boeder PDF 
in Eq. (4) and when a is large, the area under the curve of the PDF is roughly the product 
9maxP{9max) , the distribution function being sharply peaked around 9max , leading to a 
PDF approximately triangular in shape with a height P{9max) and a base of 29max- See, 
Figs. 2 and 3 to illustrate. The exponents controlling the asymptotic behaviour in Eqs. (13) 
are all consistent with respect to each other, and have been checked numerically up to quite 
large values of o; < 10^ . Note finally that the scaling results are valid in this high range of 
a values, but are inadequate in the range of ~ 1 < a <~ 500. 
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IV. SIMULATIONS AND COMPARISON WITH THE BOEDER EQUATION 



In this section, Brownian simulations are presented to obtain the PDF, and a comparison 
is made with the results of Section 3. The simulation is characterised by the presence of 
both Brownian and hydrodynamic forces acting on the macromolecule in the dilute solution. 
The hydrodynamic forces tend to act on the rod-like particles turning them in the shear flow 
with an average angular speed w, given, [3], for any orientation 9, by 

w = ^sin^{e) (19) 

where 7 is the constant shear rate of the linear flow. To simulate the effects of hydrody- 
namic forces in an interval At{n + between successive simulation events, we compute 
^(^hydiP' + 1; that is the hydrodynamic rotation about the centre of mass of the rod- like 
particles, using the following algorithm: 

Mhyd{n + 1; n) = - - ^sin\e)/^t{n + 1; n) (20) 

^QhydiP' + ^'iTi) is always positive, following the direction of the hydrodynamic flow. 
In contrast the Brownian forces in the bulk create a diffusive rotational motion of the 
macromolecules over the orientations 9. The Brownian variable A^Br(^ + l; ^) ma-Y be given 
from a random Gaussian distribution. It turns out, however, that it is more convenient to 
work with a simplified Monte Carlo random variable 

MBr{n+l]n)^±MBr (21) 

where A^sr is a fixed value for a group of simulation runs, for a given a. An appropriate 
random number generator is used in the algorithm to select random events, avoiding cu- 
mulative errors. The fixed value A^^^r must satisfy the Brownian condition corresponding 
to a simulation interval Ai(n -|- l;n) ^ l/(^rot where 5rot — ^bT /[DBrIcm\ and I cm is the 
moment of inertia of the rod-like particles about their centre of mass. Using this procedure 
/S.t{n -\- l]n) may be given as a virtual variable 

/^9l^{n+l]n) = 2DBrAt{n + l]n) (22) 
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where Dsr is the coefficient of rotational diffusion, [10]. At{n + l;n) may be ehmi- 
nated now from the algorithm. This procedure avoids ambiguities that may arise owing to 
the arbitrary choice of simulation time intervals, [5], and Eq. (20) is used to rewrite the 
hydrodynamic algorithm component as 



Typically AObt is of the order of 0.005 radians, which engenders reasonable elementary 
hydrodynamic rotations in the interval [0, vr]. This is adequate to construct viable statistics, 
[11], for the diffusion of macromolecules in solution. In Figs. 1-4, the simulation results for 
the PDF are also presented for values of a— 100, 1000, 2000, and 5000, respectively, here as 
the black dots. For any given a, the simulation results are obtained as an average over several 
simulations, each based upon 10^ hydrodynamic and stochastic events. The black dots are 
then mean values, and the error bars are the maximum uncertainty per mean value. The 
PDF simulations results arc presented in comparison with the Boeder differential equation 
solutions of the PDF, given as the continuous curves in these figures. Both are normalised to 
unity. Other comparisons of this kind between the simulations and the analytical approach 
have been made for other a values in a wide range, though not presented here. 

V. CONCLUSIONS 

The Boeder probability density function describing the average orientations of macro- 
molecules in a flowing liquid is accurately determined over a wide range of a = the 
ratio of the hydrodynamic shear rate to the rotational diffusion coefficient. The differing 
nature of macromolecules in their solutions leads to differing diffusion coefficients Dbv, 
whereas the hydrodynamic flow of the solution may be characterised by differing values 
of 7. Special analytical as well as numerical methods are developed and presented in 
order to calculate accurately this probability distribution function for a wide range of 
10"'^ < a < 10® . Scaling results are also presented valid for a > 10^ . The mathematical 
nature of the ordinary differential equation that the PDF should satisfy is revealed as a 
singular perturbation problem when a becomes larger than about 10^. 
Brownian simulations are also presented and may be shown to agree with the numerical 
solutions of the Boeder differential equation in the bulk of the flowing liquid for arbitrary 




(23) 
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values of a. This comparison confirms the sohdity of the analytical solutions. It serves 
equally as a necessary limiting reference for this type of simulation near impenetrable solid 
surfaces where the simulations become a particularly useful tool, the Boeder problem being 
analytically intractable in this region. 
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Figure Captions 

Fig. 1: The normalised probability distribution function, PDF, for the orientations of the 
macromolecular rod-like particles in a hydrodynamic flow, is presented for a high 
value of the shear to rotational diffusion ratio, a — 7^=100. The numerical solution 
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of the Boeder differential equation (continuous curve) is compared to the Monte Carlo 
simulation {black dots with error bars), for the same a 

Fig. 2: The normalised PDF, for the orientations of the macromolecular rod-like particles in 
a hydrodynamic flow, is presented for a high value of the shear to rotational diffusion 
ratio, a — ;^=1000. The numerical solution of the Boeder differential equation 
(continuous curve) is compared to the Monte Carlo simulation {black dots with error 
bars), for the same a. 

Fig. 3: The normalised PDF for the orientations of the macromolecular rod-like particles in 
a hydrodynamic flow, is presented for a high value of the shear to rotational diffusion 
ratio, a = -^^=2000. The numerical solution of the Boeder differential equation 
(continuous curve) is compared to the Monte Carlo simulation {black dots with error 
bars), for the same a. 

Fig. 4: The normalised PDF for the orientations of the macromolecular rod-like particles in 
a hydrodynamic flow, is presented for a high value of the shear to rotational diffusion 
ratio, a — ;^=5000. The numerical solution of the Boeder differential equation 
(continuous curve) is compared to the Monte Carlo simulation {black dots with error 
bars), for the same a. 
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